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We present here a formulation for the calculation of the configuration averaged lattice thermal 
conductivity in random alloys. Our formulation is based on the augmented-space theorem, intro- 
duced by one of us [1], combined with a generalized diagrammatic technique. The diagrammatic 
approach simplifies the problem of including effects of disorder corrections to a great extent. The 
approach allows us to obtain an expression for the effective heat current in case of disordered alloys, 
which in turn is used in a Kubo-Greenwood type formula for the thermal conductivity. We show 
that disorder scattering renormalizes the phonon propagators as well as the heat currents. The 
corrections to the current terms have been shown to be related to the self-energy of the propagators. 
We also study the effect of vertex corrections in a simplified ladder diagram approximation. A 
mode dependent diffusivity D~, and then a total thermal diffusivity averaged over different modes 
are defined. Schemes for implementing the said formalism are discussed. A few initial numerical 
results on the frequency and temperature dependence of lattice thermal conductivity are presented 
for NiPd alloy and are also compared with experiment. We also display numerical results on the 
frequency dependence of thermal diffusivity averaged over modes. 
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I. INTRODUCTION 

The problem of phonon excitations have been stud- 
ied extensively both theoretically and experimentally in 
mixed crystals [2, 3] as well as in insulators [4]. However, 
there are far fewer studies of lattice thermal conductivity 
in disordered alloys. Detailed comparison between theory 
and experiment on the basis of realistic models has not 
been extensive. Model calculations are mostly based on 
mass disorder, whereas in phonon problems essential off- 
diagonal disorder in the force constants cannot be dealt 
within single site mean-field approximations. Such dis- 
order cannot be ignored in a realistic calculation. 

In the past few decades, there has been considerable 
attention directed towards the theoretical understanding 
of the lattice thermal conductivity of metals. These are 
mainly due to the efforts of Klemens [5] , Ziman [6] , Call- 
away [7], Parrot [8] and others. However majority of 
these were based on model calculations, either for per- 
fect crystals or ordered alloys. Flicker and Leath [9] 
first proposed the calculation of lattice thermal conduc- 
tivity within a single-site coherent potential approxima- 
tion (CPA) using the appropriate Kubo formula. The 
single-site CPA is a mean field approximation, capable of 
dealing with mass disorder alone and is not adequate for 
treating intrinsic off-diagonal disorder arising out of the 
force constants. This was evidenced in the inability of the 
single site CPA to explain experimental life time data on 
NiPt [10]. Translationally invariant, multiple site, multi- 
ple scattering theories based on the augmented space for- 
malism [1] have recently been proposed by Ghosh et.al. 
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[11] as well as by us [12] to describe phonons in a series 
of random alloy systems : NiPt, NiPd and NiCr. These 
formalisms explicitly capture the effects of both the di- 
agonal and off-diagonal disorder. 

In this paper, we shall first introduce a Kubo- 
Greenwood type formula which relates the thermal con- 
ductivity to the (heat) current-current correlation func- 
tion. The ideas used here are very similar to those pro- 
posed by Allen and Feldman [13] except that the present 
formulation is done keeping in mind the application to a 
substitutionally disordered crystal, rather than an amor- 
phous system. For disordered alloys, configuration av- 
eraging over various random atomic arrangements have 
been carried out using the augmented space formalism 
(ASF) introduced by us [1]. The ASF goes beyond the 
usual mean-field approaches and takes into account con- 
figuration fluctuations over a large local environment. 
We shall combine the augmented space representation for 
phonons [12] with a scattering diagrammatic technique to 
get an effective heat current. This effective current con- 
sists in addition to the averaged current term, also the 
terms arising out of the disorder scattering corrections. 
We will show that these disorder induced corrections to 
the averaged current terms are directly related either to 
the disorder scattering induced self-energy matrix in the 
propagator or to vertex corrections. As far as the ver- 
tex corrections are concerned, Leath [14] had obtained 
these corrections within the framework of CPA by us- 
ing diagram summations. In this paper we shall derive 
the contribution of these corrections in a more general- 
ized context with the inclusion of diagonal as well as the 
intrinsic off-diagonal disorder arising out of the dynam- 
ical matrix. Since in an earlier communication [15] we 
have already shown that the self energy matrix and the 
Green matrix can be calculated for realistic binary alloys 
within an augmented space block recursion (ASBR) tech- 
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nique, so the present formulation will form the basis of a 
subsequent calculation of lattice thermal conductivity in 
realistic alloys. 

The rest of the paper is organized as follows. In Sec. 
II, we describe the basic tools used to calculate the lat- 
tice thermal conductivity for a crystal. In Sec. III(A), we 
briefly introduce the augmented space representation for 
phonons. In Sec. III(B,C), we derive expressions for im- 
portant physical quantities such as effective heat current, 
averaged lattice thermal conductivity and thermal diffu- 
sivity in terms of configuration averaged Green matrix 
and self energy matrix of the system. In Sec. III(D), we 
describe vertex corrections arising out of the correlated 
propagation. Sec. IV(A,B) is devoted for a description 
of the schemes for implementing this formalism to real- 
istic alloys. Few of the initial numerical results on the 
lattice thermal conductivity and thermal diffusivity for 
NiPd alloy are shown in Sec. V and an effort has been 
put forward to compare them with the available experi- 
mental data. Concluding remarks appear in Sec. VI. 



= 2v^^EE<(k)v.i^-(k)e:;,(k), 
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here 7, 7' label the various modes of vibration. 7, cjk 7' 
are their frequencies, e!^{k ),ey(k ) are the polarization 
vectors and D'^'^{'k ) is the Fourier transform of mass 
scaled dynamical matrix. 

We shall consider the case where the temperature gra- 
dient is uniform within the system. The Kubo formula 
then relates the linear heat current response to the tem- 
perature gradient field 



II. THERMAL CONDUCTIVITY 

In an earlier work [16] we had reported a formulation 
for the study of optical conductivity in disordered al- 
loys. Here wc attempt to modify that formulation for 
lattice thermal conductivity studies also in disordered al- 
loys. The formulations have similar overall similarities, 
but differ in specific details, which we would like to fo- 
cus on in this communication. The Kubo formula which 
relates the optical conductivity to a current-current cor- 
relation function is well established. The Hamiltonian 
contains a term Yliii ' A(ri,t) which drives the electri- 
cal current. For thermal conductivity we do not have a 
similar term in the Hamiltonian which drives a heat cur- 
rent. The derivation of a Kubo formula in this situation 
requires an additional statistical hypothesis [13], which 
states that a system in steady state has a space depen- 
dent Jocai temperature T(rj) = [«;B/3(rj)]~^. The ex- 
pression for the heat current has been discussed in great 
detail by Hardy [17] and Allen and Fcldman [13]. The 
readers are refereed to these papers for the details of cal- 
culation. The matrix element of the heat current in the 
basis of the eigenfunctions of the Hamiltonian is given 
by: 



s{;^,(k) = ^ (c^k^+'^ky) v^;^,(k), 



(1) 
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-00 



where 



^''-(r) = e(T) ^1^^ dX{S^{-ihX),S-'{T)), 

G(r) is the Heaviside step function, and 
S{-ihX) = e^" S e-^". 

( ) on the right-hand side of the above equation denotes 
thermal averaging over states in the absence of the tem- 
perature gradient. The above equation can be rewritten 
in the form of a Kubo-Greenwood expression 
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where, the phonon group velocity v-y.y/(k ) is given by 



E ^^"'^ qiJ. (-1, \ (If- 
d{hw^ 7) ^ 
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where (rik -y) = [e^^'^'^ -> — is the cquihbrium Bose 
Einstein distribution function and T is the absolute tem- 
perature. 

The first expression is for inter-band transitions, while 
the second expression is for intra-band transitions. For 
an isotropic response, we can rewrite the first expression 
as 



/(a;+) = lim f{w + i5), f{oj-) = lim f{w - i6). 

We have used the herglotz analytic property [12] of the 
Green operator 



^i 77' 
S^;. (k , T)5{uj' - + - '^k 7), 



where 



Si;y(k,T) = 



?i('^k 7 — Wk 7' 

We may rewrite the above equation as 



si;y(k). 



G(a; + i5) = 5ie [G(a;)] - i sgn(^) [G(a;)] . 

For disordered materials, we shall be interested in ob- 
taining the configuration averaged response functions. 
This will require the configuration averaging of quantities 
like k{zi,Z2) 

We should note that the expressions (5) and (6) are 
similar to the corresponding equations for optical conduc- 
tivity, with the heat current replacing electrical current. 
To deal with random alloys we shall follow an analogous 
procedure for optical conductivity proposed earlier [16] 



III. CONFIGURATION AVERAGING 



9m{G(k ,0;')} S''(k ,T) 9m{G(k ,uj'+lj)} 



The operator G(cj) is the phonon Green operator 
(Mw^I — $)~^. The Trace is invariant in different rep- 
resentations. For crystalline systems, usually the Bloch 
basis {|k, 7)} is used. For disordered systems, prior to 
configuration averaging, it is more convenient to use the 
basis { |k , a) }, where k is the reciprocal vector and a rep- 
resents the coordinate axes directions. We can transform 
from the mode basis to the coordinate basis by using the 
transformation matrices T^„(k ) = e"(k ). For example 



S^^(k,T) = T-^^(k)S;;^,(k,T)T;/^(k). 



If we define 



«(^i'^2) = / Tr 



S G(k ,zi) S G(k ,Z2) . (5) 



then the above equation becomes, 

-K^^{oj'-,u;'- (6) 

where 



A. The augmented space formalism for phonons 

The augmented space formalism (ASF) has been de- 
scribed in detail in several earlier papers (see [18]- [20]). 
We shall, for the sake of completeness, describe only those 
features which will be necessary for the implementation 
of our ideas in this communication. It is also important 
to introduce the notations used subsequently. The the- 
ory of phonons consists of solving a secular equation of 
the form 

(Mw^-D) u{R,w) =0, 

where Ua (R, w) is the Fourier transform of {R, t) , the 
displacement of an atom from its equilibrium position R 
on the lattice, in the direction a at time t. M is the mass 
operator, diagonal in real-space and D is the dynamical 
matrix operator whose tight-binding representations are 



M = 



D 



rriR 5a(3 Pr, 



E <n' 




Pr 



E E <n' Tnw 

R R'^R 



(7) 



where the sum rule has been incorporated in the first 
term of the equation involving D. 

Here Pr is the projection operator \R){R\ and Trr' is 
the transfer operator \R){R'\ in the Hilbcrt space H 
spanned by the tight-binding basis {j-R)}. R,R' specify 
the lattice sites and a,/3 the Gartesian directions. rriR is 
the mass of an atom occupying the position R and 
is the force constant tensor. 
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Wc shall be interested in calculating the displacement- 
displacement Green matrix G{R,R',w'^) 



G{R,R',w^) = {R\ (Mw^ 



D 



\R'). 



Let us now consider a binary alloy A^By consisting 
of two kinds of atoms A and B of masses tha and 
niB randomly occupying each lattice sites. We wish 
to calculate the configuration-averaged Green matrix 
<C G{R,R',w'^) :s>- We shall use the augmented space 
formalism to do so indicating the main operational re- 
sults here. For further details we refer the reader to the 
monograph [19]. The first operation is to represent the 
random parts of the secular equation in terms of a ran- 
dom set of local variables {n^} which are 1 if the site R is 
occupied by an A atom and if it is occupied by B. The 
probability densities of these variables may be written as 



Prinn) = x S(nji-l) + y S{nji) 

= (-I/tt) 9m(Tfl| (riRl-Nfi)-' \U}, (8) 

where x and y are the concentrations of the constituents 
A and B with x + y = 1. Nr is an operator defined on 
the configuration-space (pu of the variable ur. This is of 
rank 2 and is spanned by the states {ITij), Ufl)}i 



M = A(m) 7(g) 7-1- B(m) ^p^OPfl 

R 

+F(m) J2 Tr ® Pr 

R 

= <CM> +M', (10) 



where 



A(X) X > = (a;X^ + yXg), 
B(X) = iy-x) (Xa - Xb), 
F(X) = ^{Xa- Xb). 

Similarly the random off-diagonal force constants ^^/J' 
between the sites R and 7?' can be written as 



= <A«HnH'+<s(l-««)(l-riii') + 

nR{l -nR')+ nR'{l - ur) 



^BB 



a/3 



a/3 \ 



BB-^'^AB nRnR,+ 



(<s-<s) (nR + riR,). 



(11) 



U) = V^\Ar)+^\Br), I Ir) = ^\Ar)-V^\Br), Let us define the following : 



Nr 



xPr 



VPr- 



Ti 



at 

~R ' 



_ ^a/3 , ^a/3 _ cyff^o 
^(2) ~ ^AA^^BB ^^AB- 

In augmented space the off-diagonal force constant ma- 
trix becomes an operator 



where \Ar) is the state in which an atom of the type 
A occupies a site 7?. \Br) is similarly defined, pjj = 
I T-r)(T-R I and = | | are projection operators 



and tH = I ]R){iR I and = \ ifl)(Tij | are transfer 
operators in the configuration space 4'r. 

In terms of random variables ur, the mass operator 
can be written as 



M = ^ + riR {6m) 5ap Pr, 

R 



6m = mA—ms 
(9) 

According to the augmented space theorem, in order to 

obtain the configuration-average wc simply replace the 
random variables ur by the corresponding operators Nr 
associated with its probability density, as in Eqn. (8), 
and take the matrix clement of the resulting operator 
between the reference states. For a full mathematical 
proof the reader is referred to [19]. 



Ur — >Nr = xi + {y - x) p^ji + ^/xy T^K 
Using the above we get, 



RR' 

^{y - X) (p^T^/ + p^^,Tl^) + xyTl^Tl} } 

= E « ^RR' » ~^ ® + E ® • 



RR' 



RR' 



RR' 



The sum rule gives the diagonal element. 



wL. = -E^ E ^^tR'»i)®PR 



(dia) 



R I R'^R 



E \i:^fn')®PR 

R I R'^/tR 
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The total dynamical matrix in the augmented space is 



D = «D»-^ ^ ^1% 

R yR'i^R 

+ E "^"RR'^TnR, 

RR' 

= <C D » + D'. 



Pr 



(12) 



where 



(14) 



Referring back to Equations (10), (12) and (13) we get 



<G(i?,i?',w;")» 



we define 



R [ R'^R 
R R'^R 



'Pr + 



with 



Ti 



B(m) p\j^ + F(m) 

^'rI {p'r+p'r')+^Tr' + 
^'Zp'rp'r^+^Z {p'rT^^+T^'p'r-) 
+ (16) 



The boldface operators are 3x3 matrix representations where 
in the three Cartesian directions. 

The augmented space theorem [1] now states 
that the configuration-average of the Green matrix 
<C G{R,R' ^vj^) may be written as 

^G{R,R',v?) > 

= ({0}0i?| (mm;2-d)"' \{%}(S)R'), (13) 

where M and D are the operators which are constructed 
out of M and D by replacing all the random variables 
ur (or nn' ) by the corresponding operators (or Nn' ) 
as given by Eqn.(lO) and (12). These are the opera- 
tors in the augmented space Vl = 7i $5 The state 
\R® {0}) is a state in the augmented space, which is the 
direct product of the real-space and the configuration- 
space bases. The configuration-space $ = YI'r. 4'R. is of 
rank 2''^ for a system of N-lattice sites with binary distri- 
bution. A basis in this space is denoted by the cardinal- 
ity sequence {C} — R2, ■ ■ ■ , Rc\ which gives us the 
positions where we have a 1 1) configuration. The config- 
uration {0} refers to a null cardinality sequence i.e. one 
in which we have | t) at all sites. 

The virtual crystal (VGA) Green matrix is 
g{R,R',w'') = ({0}®i?|(<M» ^g^-)-i||0|^^/^^ 



D(i) 
D(2) 

D(3) 
D(4) 

d(5) 



(y-x) $(1), 

(y-x)2 $(2), 

Vxy (y-x) $(2), 
xy *(2)- 



Scattering diagrams are obtained by expanding the Eq. 
(15) as an infinite series and the terms B, F and D^^-* to 
D^^^ are represented as scattering vertices (see Fig. 1). 
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J3) 



^ ■ W 



+ M'a;2-D')"' |{0}®J?') 
m^R\ (g-i-Di)"'|{0}®i?'), 



J5) 



FIG. 1: The scattering vertices for the averaged Green func- 
tion 

(15) 

It might be instructive to understand what these scat- 
tering vertices represent physically. If we look at Eq. 
(15), we note that the term Di leads to the creation or 
annihilation of conGguration fluctuations over and above 
the virtual crystal description. The vertices F shown in 
Fig 1, create and annihilate a configuration fluctuation 
at a given site because of mass disorder, while the ver- 
tex B counts the number of such fluctuations at a given 



6 



site. These are the only type of configuration fluctu- 
ations we can have if we had single site mass disorder 
alone. The single-site mean-field approximations like the 
single-site coherent potential approximation (ICPA) can 
ideally deal with situations where we ignore the other 
vertices in Fig 1. 

The vertices D^^^ also describe creation and annihila- 
tion of configuration fluctuations at single sites. That 
is, fluctuations at any one end of the two-site dynamical 
matrix. Similarly, the vertices D^^) count the number of 
configuration fluctuations at any one of the two ends of 
the dynamical matrix. These arc also single-site config- 
uration fluctuations but arise due to fluctuations in the 
two-site dynamical matrix. These may also be treated 
with some variant of the ICPA. For example, there are 
versions of the ICPA which assume 2D^^ = + 
D^^. With such an assumption only the single site con- 
figurations fiuctuation vertices are non-zero. 

The vertices D*^^) describe creation and annihilation of 
two configuration fluctuations, one at either end of the 
two-site dynamical matrix. The vertex D^''^ counts the 
number of configuration fluctuations at either end of the 
dynamical matrix. The vertices D^^' are mixed types 
which both create or annihilate a configuration fiuctua- 
tion at one end of the dynamical matrix and count the 
number of fluctuations at the other end. 

These last three vertices describe configurations fluctu- 
ations which are essentially two-site and cannot be prop- 
erly described within a single-site mean-field approxima- 
tion. 



Scattering diagrams for the averaged Green 
functions and mean-field approximations 






In our earlier paper [15] we had developed a multiple 
scattering picture for the configuration averaged Green 
function. The idea is very similar to that of Edwards 
and Langer [21, 22] in the context of purely diagonal dis- 
order. The scattering vertices associated with the terms 
in the "perturbing" part of the dynamical matrix : Di 
are shown in Fig. 1. 

The end point of that formulation was the derivation 
of a Dyson equation : 



« G >= g-Kg s «: G > . 

For homogeneous disorder we have shown earlier that we 

have translational symmetry in the full augmented space 
[23]. We can then take Fourier transform of the above 
equation to get 



« G(k , E) >= g(k , ^)+g(k , E) S(k , E) < G(k , E) > 

The diagrams for the self-energy are skeleton 
diagrams [34] which have the structure as shown in Fig. 
2. Each of the 25 different diagram starts with any one 
of F, D^^^ or D*^^^ and the central dark semicircle repre- 
sents all possible arrangements of scattering vertices to 
all orders. 

Let us now examine some specific Edwards-Langer 
scattering diagrams, in some detail, in order to under- 
stand their physical significance and relation to mean- 
field approximations. The first three diagrams on the 
first row of Fig 3 arise because of purely diagonal dis- 
order in mass. Of these the first two diagrams describe 
self-energy corrections to the VCA propagator because 
of configurations fiuctuations at a single site. These dia- 
grams are closely related to the diagrammatic treatment 
discussed by Leath [24]. Referring to that earlier work, 
we note that these diagrams are explicitly included in the 
ICPA. The self-energy arising out of such diagrams is 
site-diagonal, or k-independent in reciprocal space. The 



R RRRRRR'RR' 



R'R R RR' RR' R R'R 



R R RR' R'R R R'R RR" 



R R' R R' RR' R R'R" R R' R R' RR' R'R 



FIG. 2: Structure of the skeleton diagrams for the self-energy 



FIG. 3: Details of some skeleton diagrams for the self-energy 
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third diagram is the smallest order diagram of this type 
which describes joint configuration fluctuations of two 
distinct sites. These type of diagrams take us beyond 
the ICPA. For the ICPA we ignore these fourth order 
diagrams, and all higher order diagrams, to all orders, 
describing joint configuration fluctuations of more than 
two sites. Within this approximation, the first inaccu- 
rate moment of the density of states is of order eight. 
If we include diagrams to all orders which describe joint 
configuration fiuctuations of two sites we are lead to the 
2CPA. This has been described in detail in the work of 
Aiyer et.al. [25] and Nickel and Krummhansl [26]. 

The diagrams in the second row of Fig 3 describe self- 
energy corrections due to configuration fluctuations at 
a single site for both diagonal (mass) and off-diagonal 
(force-constant) disorders. The first and the third dia- 
grams lead to a diagonal self-energy, while the second 
and fourth diagrams lead to an off-diagonal contribution 
to the self-energy in real space. For off-diagonal disorder 
even these second order diagrams can lead to off-diagonal 
(real space) or k-dependent (reciprocal space) self-energy. 
Ignoring these contributions ( to all orders) will lead to 
a ICPA type approximation where even the fourth mo- 
ment of the density of states will be inaccurate. This 
had been noted earlier for ICPA in off-diagonal disorder 
problems. The diagrams in row three of Fig 3 are very 
similar contributions from configuration fiuctuations at a 
single site but arising out of pure off-diagonal disorder. 

The diagrams on the last row of Fig 3 describe self- 
energy contributions coming from joint fluctuations of 
two sites arising out of off-diagonal disorder. All such di- 
agrams take us beyond the ICPA and some of them con- 
tribute to a self-energy which is off-diagonal (real space) 
or k-dependent (reciprocal space). 

The formal summing up, to infinite order, of diagrams 
which involve configuration fiuctuations involving single 
sites has been discussed by Leath [24]. This leads to the 
ICPA and within this approximation the first inaccuracy 
occurs in the eighth moment of the averaged density of 
states. Nickel and Krummhansl [26] have discussed in 
detail how to sum up, to infinite order, diagrams which 
also include joint configuration fluctuations of two sites. 
It is clear from this paper that any generalization of a 
purely diagrammatic treatment of joint multi-site config- 
uration fluctuations is a very difficult task indeed. How- 
ever, a recursion based calculation of the self-energy, as 
proposed by us in our earlier paper [15], includes the con- 
tribution of such diagrams. We note that if we calculate 
the continued fraction coefficients exactly up to n steps 
then all the first 2n moments of the density of states are 
exact. In addition the Beer-Pettifor terminator ensures 
both the herglotz analytic properties of the Green func- 
tion and accuracy of the asymptotic moments. As in our 
earlier papers [12, 15], therefore, we propose the aug- 
mented space recursion as an alternative method for the 
calculation of averaged Green functions including effects 
of joint multi-site configuration fluctuations. 



C. Generalized scattering diagrams for the 
averaged thermal conductivity 

We now go back to equation (5) and discuss the con- 
figuration averaging of the two-particle Green functions 
of the kind k{zi, z-i). The augmented space theorem im- 
mediately implies that 



< k(zi, Z2) ^ 



87r3 



Tr 



SG(k ,^1) SG(k ,02) {0} 



In real space the random expression S for a binary 
alloy can take the values S^^-^, S^^, S^-^ or S^^ with 
probabilities , xy, yx and respectively. We may 
rewrite the current Siioi,R'i3 as 



+ ^RiR'0 (1 - riR) nw + ^rZr'0 (1 - - "flO- 

Following the same procedure as for the single particle 
Green functions we get 



S = 



EE 

Ra R'a' 



(17) 



where 

s(i) = 
s(^) = 

and 

S(i) 

^Ra,R'l3 
S(2) 



(t/-x)S«, s(2) = 
{y-xfS^^\ S(4) = 
xy S(2). 



^ ' ^Ra,R'(3 ~ °Ra,R 



xy SW, 

xy{y-x) §(2), 



'0) - V ( 



cBB 

^Ra,R'l3' 



^BA \ 
^Ra,R'l3 J 



cv^; aAA I cBB Sab aBA 

^Ra,R'l) ~ ^Ra,R'l3 '^Ra,R'l3 ^Ra,R'0 '^Ra,R'l3- 

Equation (17) is very similar to (16), which shows the 

terms arising in augmented space due to the disorder in 
the dynamical matrix. Fig. 4 shows the sixteen different 
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<s> 



/ \ ^ . 

\/i V/ / 



FIG. 4: The scattering vertices associated with the random 
current terms 



scattering vertices arising out of Eq. (17). We may com- 
pare the scattering vertices arising out of disorder in the 
current term with those that arise from disorder in the 
dynamical matrix shown in Fig. 1. We note the close sim- 
ilarity between them leading to both terms which arise 
because of configuration fluctuations at a single site and 
also two sites. 

Note that the averaging of k{zi,Z2) involves averag- 
ing of a product of four terms two involving the thermal 
currents S and two Green functions. These random func- 
tions are not independent, but are all intrinsic functions 
of the random occupation variables nij. The average of 
the product is therefore not the product of the averages. 




FIG. 5: Different classes of scattering diagrams for the ther- 
mal conductivity 



Let us now examine the various classes of scattering 
diagrams. The simplest one is the type shown on the 
top left of Fig 5. These diagrams involve two averaged 
current terms and two propagators decorated with scat- 
tering diagrams in all possible ways. The contribution of 
such a term (the zero-th order term) is then, 



Tr 



8^3 

< S » G(k , Z2) » 



< S G(k ,zi) > 



The next type of scattering diagrams are the types 
shown on the top right of Fig 5. These diagrams connect 
a current term to a propagator. The first type of scat- 
tering diagrams are those in which disorder lines do not 
connect the two phonon propagators. Fig 6 shows the 
general structure of such diagrams. It is clear from the 
figure that these diagrams renormalize the current terms. 
Fig 7 shows all possible renormalized currents arising out 
of diagrams of this class. 






FIG. 6: Two examples of scattering diagrams where no dis- 
order line joins the two phonon propagators 



In our earlier paper on optical conductivity [16] we had 
shown that this type of diagrams provide the predomi- 
nant disorder correction to the current terms in the elec- 
tronic problem. If we compare them with the diagrams 
for the self energy (Fig. 2) we note that for the diagrams 
in the left panel of Fig. 7, the only difference is that the 
leftmost scattering vertex is replaced by a very similar 
current term. Similarly the diagrams on the right panel 
are the same as the self-energy diagrams except that the 
last vertex is replaced by a very similar current term. In 
all these diagrams of Fig. 7 the left and rightmost diag- 
onal term similar to the vertex F of Fig. 1 is of course 
missing. We may then go on to write 
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3i 






FIG. 7: Scattering diagrams contributing to effective heat current 



ASi(k)(^i,^2) = 

2 (S(2)(k ) + S(5)(k )) [A(k , z,)] S(k , zi) . 
S(k , ^2) [A(k , ^2)] 2 (S(2) (k ) + S(^) (k )) 

where 



(18) 



A(k , z) = + 2D(2)(k ) + 2D(5)(k ). 

The next most dominant disorder corrections come 
from the diagrams shown in the middle row of Fig 5. A 
general such diagram is shown in Fig. 8. These diagram 
connect one current term to two propagators. These dia- 
grams also renormalize this current term. Again looking 
at these diagrams we note that they are also related to a 
self-energy diagrams with vertices at both ends replaced 
by currents. We therefore may write 



AS2(k )(zi,Z2) = 
S(k ,Z2) fA(k ,Z2) 



S(5)(k)[A(k,zi)l S(k,zi) 



The effective current is then given by 



Scff(k ,2:1,22) = 
< S(k ) > +ASi(k ){zuZ2) + AS2(k ){zi, Z2). (20) 



In our earlier paper [16] on a similar electronic prob- 
lem we have argued that these are the dominant disorder 
corrections to the current term. With these corrections 
we obtain 



K(1){Z1,Z2) 



8^ 



Tr 



Seff(k ,Zi,Z2) < G(k,Zi) » 



SlfiOsi ,zi,Z2) < G(k, Z2) > 



(21) 



+S(k,2i) [A(k,2i)l S(5)(k) [A(k,22)l S(k,22 



FIG. 8: 
current 



Since averaging the thermal conductivity involves av- 
eraging the product of two propagators and two current 
(19) terms, we may physically understand the different terms 
in Fig. 5 in the following way : the top left diagram is 
related to the product of the averages. The term on the 
top row to the right is related to a pair of joint averages 
of one propagator and one current term. The diagrams 
on the second row are related to the product of the aver- 
age of one current and the average of the product of two 
propagators and a current. The bottom most diagram is 
More scattering diagrams contributing to effective related to the product of two averaged currents and the 

average of the product of two propagators. 
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D. Averaged thermal diffusivity 

For a harmonic solid, a temperature independent mode 
diffusivity is defined as 

(k ) = TT ^ -^Si;^'(k )S!;,^(k ) (5(a;k , - y). 

This is an intrinsic property of the 7-th normal mode and 
provides an unambiguous criterion for localization. 

The averaged thermal diffusivity (averaged over 
modes) is then given by 



7 



8^ 



(22) 



Assuming isotropy of the response, we can rewrite the 
numerator of above equation as 



7 7' 



S{u)' — Wk 7')'^('^k 7 — U)')5{u) — Wk 7), 



where 



Wk 7 



s^;y(k ). 



We may again rewrite the above equation for Dt^t as 



fk 

87r3 



Tr 



Sm{G(k ,w')}S''(k 



9TO{G(k ,w')}S^(k )3?m{G(k ,u;)} 



The averaged thermal diffusivity can then be expressed 
as ( for an isotropic response ) 



D(a;) = l^D^^iu;) 



TT 

sTcFk 



^Tr[9m{G(k,u;)}] 



(23) 



For disordered material, we shall be interested as before 
in obtaining the configuration averaged thermal diffusiv- 
ity. We have already discussed the configuration aver- 
aging of the two particle Green function using scatter- 
ing diagram technique in Sec. III(C). It has been found 



that the net effect is to replace the current terms by an 
effective heat current Se//(k ,zi,Z2)- The effective cur- 
rent is a sum of average current and the terms arising 
out of the disorder correction. As in the case of optical 
conductivity calculation [16], it will be shown that the 
overall contribution of disorder correction terms to the 
thermal conductivity is very small as compared to the 
average current <C S(k ) Keeping in mind this effect 
of disorder correction terms to the heat current, the con- 
figuration averaged thermal diffusivity can be expressed 
(to a 1st order approximation) in the form 



<C D(a;) >= - 



E„ « nzi^) » 



/ -^Tr[3m< G(k ,0;) »] 



(24) 



where 



«A';'^H»^ ^[dio'f ^Tr 



9m < G(k , w') » 



«; S''(k ) {G(k S''(k ) Sjm{G(k ,0;)} > 



(25) 



E. The vertex correction 

We shall now examine the scattering diagrams we have 
left out, namely, those in which disorder lines connect 
both the propagators directly. One such diagram is the 
bottom-most diagram shown in Fig 5. These lead to 
vertex corrections due to correlated propagation. In gen- 
eral we obtain a Bethe-Salpeter equation for the averaged 
two-particle propagator. This is diagrammatically shown 
in Fig 9. 





FIG. 9: The Bethe-Salpeter equation for the thermal conduc- 
tivity 

To obtain the vertex correction, shown in Fig 9, we 
shall consider only a special class of vertex diagrams in 
this communication : namely the scattering diagrams 
which arc built out of repeated vertices. These are called 
the ladder diagrams and can be summed up to all or- 
ders. This is the disorder scattering version of the ran- 
dom phase approximation (RPA). 

A few amongst all the possible scattering diagrams for 
the ladder type of vertex corrections involving the ver- 
tices B, F, D(i), D(2), d(3), d(4), D(5) are as shown in 
Fig. 10 
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FIG. 10: The scattering diagrams for the ladder type of vertex 
correction. 



The contribution of diagrams from the first hne of Fig. 
10 in the four legged vertex (shown in the extreme left 
column of the same figure) is given by 



1st line 



+2zt i)g F« + AD':> D^^>. (26) 



(2) r.(2) 



The contribution of the second line of diagrams and all 
those obtained from them by simple reflection operations 
about the vertical and horizontal directions is given by 



2nd line 



(-Pai/"5ai/' ) Giii,',Rvlt {Fv'i fjSv" jj) 



(-Fj/i/'^i/i/') Grv'^Rvi, {F,y»sSv"s) 



jj' v" 



i?i'j(^2'B.5<5.5)l + 2 [(^i2s„^<5„^)Z?|,yl + {z^B^ii8^fi){z2^B,M 



+ 2 '^^^{zx'Fav'^av') Giiv',Rvii D^^,] 



+ 2 

+ 4 



J2 Gr.',r.„ {zi^F,„p6,„p) 



{Z2^B^S5^S) + 2 {zi-Bafi^al}) 
{z2^B^s5u&) +2 {zi-BapSajs) 



^(22'i^..'5..') Gr,,,r,„ dI^} 
v'v" 

Gr,,,r,„ iz2^F,„sS,"6) 



E Gr.,,r.„ D^, 



+ Df^^z^'B^M + {zx'B^^bocfi)D'°l. 



{Z2 Bj^sS^s 



(5) 



(27) 



The contribution of rest of the diagrams may be obtained 
in exactly the same way. 

If the sum of all possible scattering diagrams contributing 
to the four logged vertex shown in the left column of Fig. 
10 is denoted by W^^^ then the infinite scries of ladder 
diagrams shown in Fig. 11 may be directly summed up 
as follows : If we define 



r(^i,22) = 






Jbz 


T{ZI,Z2) = 


j 




Jbz 


1^{Z.,Z2) = 


j 

JBZ 



87r3 
fk_ 

87r3 

87r3 



G(k,02)S<'*'(k,0i,^2)G(k,^i), 

G(k ,Zi)S<=«(k ,Zi,Z2)t G(k ,Z2), 
Ga/3(k ,Zi) Grysi}i ,Z2), 



Then 
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salient features of the technique. In the block recursion, 
we start from a matrix basis of the form where 
J is the discrete labeling of the augmented space states 
FIG. 11: The structure of infinite series of ladder diagrams. and the a, (3 labels Cartesian directions. 

For a reciprocal space based calculation (as in the 
present case) we start with 




A{zi,Z2) = W + W0W + W0W0W+... 

= W(0i, Z2) (l - 0(^1, ^2)W(^i, ^2)) ~' . 

(28) 



Thus we get 



> 



a.0 vS 



: Tr [r(zi, 2:2) ® f{zi,Z2) © A{zi,Z2) 



(29) 



The configuration average of the two-particle green 
function is then given by 

< k{zi,Z2) >=« K(^l){Zi,Z2) > + < Ak{Zi, Z2y'"^'^'''' > 

(30) 



IV. DETAILS OF NUMERICAL 
IMPLEMENTATION 



Although we have used the scattering diagrammatic 
technique to analyze the effects of disorder scattering on 
the thermal conductivity and obtain relation between the 
effective current and self-energy, we shall not use this ap- 
proach to actually numerically obtain the thermal con- 
ductivity for a real alloy system. If we look at the earlier 
sections we note that what wc need to obtain are essen- 
tially the configuration averaged Green matrices and the 
matrix self-energies. For this we shall use the augmented 
space block-recursion [15] and also the Brillouin zone in- 
tegration scheme developed by us earlier [27]. In this 
section we shall briefly discuss the two techniques. 



*Si/3 = Uipj,, + u!^^Sj,2, 

where 

[A{m 1)] [A{m ^)\ 

with 

A{X) = {xXa + vXb), F{X) = {y- x){Xa - Xb). 

The remaining terms in the basis are recursively ob- 
tained from 



E^(n + l) „(n-Hl)t 



(n) 



J' 0' 

E 

0' 



(f,{n) An) _ ^ ^ 



V 7 -Dfl/fl, 



where H is the lattice vibrational Hamiltonian (under 
harmonic approximation) in the augmented space. We 
then use the Gram-Schmidt orthonormalization method 

to keep track of orthogonality condition among the vari- 
ous cohimns of ^^t"''^ and hence to calculate the matrix 

coefficient 3^"+^'. Since the Hamiltonian in this new 
basis is block tri-diagonal, so the green matrix can be 
expressed as a matrix continued fraction in terms of the 
coefficients A",B("+^^ as follows : 



Q(n) = I- A^") - b("+^^''' g("+^) b("+^) 



(31) 



We calculate the matrix coefficients up to a n = A^o 
and approximate at coefficients > TVq by A and B.We 
then write fov a, N ^ Nq, 



Block Recursion 



g(^) = 



(^2 - i6)I 



As is clear from the expression of configuration aver- 
aged lattice thermal conductivity and thermal diffusivity 
that the numerical evaluation of these quantities require 
the full Green matrix and the self energy matrix instead 
of their diagonal elements alone. 

In a recent communication [15], we have already de- 
scribed the block recursion technique, which calculates 

the entire green matrix and the self energy matrix. For 
the sake of completeness, we shall describe in brief the 



and then iterate 

= \w^I - A - Bt g("+i) b] 

for n > No. 

The self-energy follows from the Dyson equation 
E = - G-^ 



The Green matrix and self-energy matrices are essen- 
tial inputs for the calculation of the thermal conductivity 
and effective currents. 



B. Brillouin zone integration for disordered 
systems 



The need of efficient techniques for Brillouin zone (BZ) 
integration in solid state physics has been widely appre- 
ciated in recent years. Such techniques are of great im- 
portance in the numerical calculation of density of state, 
conductivity, susceptibility, dielectric function etc. The 
tetrahedron method allows us a very accurate k-space 
integration for both the spectral functions and conduc- 
tivities for ordered systems. Recently our group has de- 
veloped a generalization of this technique for disordered 
alloys. The spectral functions are now no longer delta 
functions, but Lorcnzians with a disorder induced non- 
zero half- widths [27]. We will use the efficient codes de- 
veloped by our group to carry out the integrations over 
the Brillouin zone. We refer the reader to the above ref- 
erenced paper for details of the calculation. 



V. RESULTS AND DISCUSSION 

We have calculated the lattice thermal conductivity of 
NiPd alloys. Fig 12 shows the thermal conductivity as a 
function of the frequency at T=170 K. We shall analyze 
the contributions of various disorder correction terms to 
the averaged thermal conductivity at this temperature 
from Fig 13. The figure shows separately the corrections 
due to disorder renormalization of the current terms and 
that due to the vertex correction. We note that the major 
contribution comes from the zero-th order term involv- 
ing only the averaged currents and averaged propagators. 
The largest correction comes from the term which arises 
due to the disorder renormalization of the current terms, 
while the vertex correction in the ladder approximation 
also has a small but non-negligible contribution. 

Direct comparison with the experimental data on these 
systems is difficult, because the experimental thermal 
conductivity also has a component arising out of the con- 
tribution from electrons. Figure 14 shows the tempera- 
ture dependence of lattice conductivity. The top panel 
shows our theoretical result for the iViggPdoi alloy at 
three different frequencies. The bottom panel shows the 
experimental data [29] on the total ( electronic + lat- 
tice ) thermal conductivity of the same 99-01 NiPd alloy. 
Since the frequency is not mentioned in the experimen- 
tal data, we assume that it must be for low frequencies. 
The best comparison then will be between the middle 
(black) curve on the top panel and that in the bottom 
one. The two agree qualitatively, except at low tem- 
peratures where we expect the electronic contribution to 
dominate. In order to understand whether the devia- 



20 - 



A 10 
A 

y 

V 
V 

5 



1 


1 1 1 




— «K» (VCA + All Collections) - 




\ ForNi Pd alloy 

\ 50 50 


/ 


\ T=170K 




1 — 1 



Frequency V(THz) 

FIG. 12: Thermal conductivity as a function of frequency for 
NisoPdso alloys at T=170 K 



tion does arise from the electronic contribution, we have 
compared the top panel with the thermal conductivity 
of amorphous-Si [28] , shown in the middle panel. In a-Si 
the electrons near the Fermi level are localized and hence 
cannot carry any current. Almost the entire contribution 
should come from the phonons. The behaviour of the two 
panels are quite similar. 

A careful inspection of our results (top panel of Fig. 
14) indicates that at low temperatures where only low- 
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FIG. 13: (Color Online) Corrections arising due to the dis- 
order renormalization of the current term and the ladder ap- 
proximation to the vertex correction 
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FIG. 14; (Color Online) Thermal conductivity vs temperature 
T(K) for NiPd alloys and Amorphous Si. The top panel shows 
our results on the lattice conductivity for NiggPdoi alloy at 
three different frequency cut-off ly. The middle panel shows 
the lattice conductivity for amorphous Si [28] at three differ- 
ent cut-off frequency, while the panel at the bottom shows the 
experimental data [29] for the total thermal conductivity {— 
lattice -I- electronic contribution) of the same NiggPdoi alloy. 



energy vibrations are excited, k(T) is approximately a 
quadratic function of T. Figure 15 shows a plot of n 
vs T^. The calculated curve fits reasonably well with a 
straight line. This has been seen in experimental obser- 
vations [29] . Additional scattering processes leading to a 
different temperature dependence of lattice conductivity 
become apparent at higher temperatures. At T >25 K, 
k(T) rises smoothly to a T-independent saturated value. 
The dominant mechanism in this regime is the intrin- 
sic harmonic diffusion of higher energy delocalized vibra- 
tions. These modes have not been well described by most 
previous theories. 

The middle panel for amorphous Si has the similar 
qualitative behaviour for k(T) as ours. The three curves 
in this panel stands for three different frequencies. It is 
clear from the two panels (top and middle) that the sat- 
uration of k{T) for NiPd alloy starts at an earlier tem- 
perature as compared to that for amorphous Si. The 
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FIG. 15: K vs for low temperatures for 50-50 NiPd alloy 



magnitude of lattice conductivity in NiPd alloy is also 
larger than the value in amorphous Si, which is obvious 
because NiPd alloy is a metallic system with larger heat 
conduction. 

In the bottom panel, which shows the experimental 
data for total thermal conductivity k = Ke + Hg-, where He 
is the electronic contribution and Kg the lattice contribu- 
tion to the thermal conductivity, one can notice various 
similarities with the results of top panel. For example 
both the panels have (i) an approximate T^-dependence 
[29] in the low temperature regime and (ii) the value 
of k(T) approaching a T-independent final or saturated 
value at higher temperature side. In spite of these simi- 
larities, the overall agreement between the top and bot- 
tom panels is not satisfactory. This is mainly due to the 
existence of a broad peak like behaviour in the experi- 
mental curve unlike the theoretical results. This feature 
has been described in detail for a series of dilute alloys 
like Ni-Co, Ni-Fe and Ni-Cu in [29]. We would like to 
ascribe this difference between the theoretical and ex- 
perimental data to the extra electronic contribution to 
thermal conductivity hidden in the experimental curve. 

The thermal diffusivities D(z/) are important because 
the effect of disorder is often manifested in them more di- 
rectly than in the conductivities. Not only that the ther- 
mal diffusivity also give an approximate idea about the 
location of mobility edge as well as the fraction of delo- 
calized states. Figure 16 shows the frequency dependence 
of diffusivity for NisoPdso alloy. There are basically two 
regions of large thermal diffusivity, one near the lower 
frequency region (~ 0.5THz) and the other around a 
somewhat higher frequency region (~ 1.25THz). Above 
3 THz, there is a smooth decrease of diffusivity approx- 
imately linear in frequency. The Fig 17 shows a plot 
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FIG. 16: The configuration averaged thermal difFusivity D{i/) 
vs frequency v for NisoPdso alloy. The inside panel shows the 
phonon density of states for the same alloy. 



of D(t/)cx (i^c — t^)"- At the critical frequency I'c— 6.3± 
0.014, D{i') vanishes to within a very small level of noise. 
The critical exponent a ~ 1.0041. A knowledge of these 
critical parameters leads to the information about the lo- 
cation of mobility edge. The critical exponent of a ~ 1 in 
our case agrees with scaling and other theories of Ander- 
son localization [30] and the critical frequency i^c locates 



Approximate Linear decrease of Diffusivity (D) above v =3 THz 
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FIG. 17: (Color Online) The behaviour of thermal diffusivity 
as a function of {vc — v)" 



the mobility edge above which the diffusivity is zero in the 
infinite size limit. Once the mobility edge is located, the 
fraction of de-localized states may be obtained by evalu- 
ating the area under the D{v) vs v curve from v — to 
V = Vc- For Ni5oPd5o alloy ( see Fig. 16 ), the mobility 
edge is located at Vc — &.'6\ATHz. The panel at the top 
right corner of Fig. 16 shows the phonon density of states 
for the same alloy. The shaded region in this panel gives 
the information about the localized states. The fraction 
of localized states is actually the area beyond the critical 
frequency Vc- 



VI. CONCLUSIONS 

We have formulated a theory for lattice thermal con- 
ductivity based on a realistic model. The augmented 
space approach allows us to go beyond the standard 
single-site mean-field theories and the augmented space 
recursion allows us to include, in the calculation of our 
averaged propagators, effects of joint fluctuations at more 
than one site. The augmented space approach has ear- 
lier been generalized to include effects of short-ranged 
ordering [31, 32] as well as local lattice relaxations due 
to large size mismatch of the constituents [33] . Calcula- 
tions including these effects have been earlier carried out 
for the electronic case. We propose to apply the same 
technique to phonons in disordered alloys. The scatter- 
ing diagram approach proves to be useful in analyzing 
and calculating the disorder corrections to the averaged 
current. These are shown to be the dominant corrections 
and are related to the self-energy. Next in importance, 
we have studied the effect of vertex corrections arising 
out of the correlated propagation. We have shown ex- 
plicitly how to obtain these corrections within the ladder 
diagram approximation. Our formalism explicitly takes 
into account fluctuations in masses, force constants and 
heat currents between different nuclei. For the calcula- 
tion of the averaged propagators themselves we have used 
the augmented space recursion with the Beer-Pettifor ter- 
minator scheme. Our efficient Brillouin zone integration 
codes for disordered alloys makes the numerical calcu- 
lation stable and accurate. We have already shown in 
an earlier communication [15] that the approximation 
involving termination of the matrix continued fraction 
expansion of the green matrix retains the essential her- 
glotz analytic properties of the diagonal green function. 
Our numerical results on the temperature dependence of 
lattice conductivity favours a general trend of other the- 
oretical results as well as the experimental data. 
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